## plot r-star estimates from UC model

rm(list=ls())
source("R/data_fns.R")
graphics.off()
library(dplyr)
q <- 0.68   ## 68 percent credibility intervals (1 SE width for normal distribution)

cat(100*q, "percent CIs\n")
cat("Lower quantile:", lq <- (1-q)/2, "\n")
cat("Upper quantile:", uq <- (1+q)/2, "\n")

load("results/uc_estimates_y1.RData")
stopifnot(all.equal(data$rstar.mean, colMeans(X[,1,])))
data <- data %>%
    rename(r1 = rr,
           rstar.r1 = rstar.mean) %>%
    mutate(rstar.r1.lb = apply(X[,1,], 2, quantile, lq),
           rstar.r1.ub = apply(X[,1,], 2, quantile, uq)) %>%
    select(year, r1, starts_with("rstar.r1"))
env <- new.env()
load("results/uc_estimates_y10.RData", envir=env)
d10 <- env$data %>%
    rename(r10 = rr,
           rstar.r10 = rstar.mean) %>%
    mutate(rstar.r10.lb = apply(env$X[,1,], 2, quantile, lq),
           rstar.r10.ub = apply(env$X[,1,], 2, quantile, uq)) %>%
    select(year, r10, starts_with("rstar.r10"))
data <- left_join(data, d10)
data %>% filter(!is.na(r1)) %>% pull(year) %>% range
data %>% filter(!is.na(r10)) %>% pull(year) %>% range

## external r-star estimates
data <- data %>%
    mutate(yyyymm = year*100+12) %>%
    left_join(loadDN()) %>%
    left_join(loadJM()) %>%
    left_join(loadLW()) %>%
    left_join(loadKiley()) %>%
    select(-yyyymm)
nms <- c("rstar.dn.sm", "rstar.jm.sm", "rstar.lw.sm", "rstar.kiley.sm")
data$rstar_mean <- rowMeans(as.matrix(data[nms]))

## Figure 1 - rstar from UC models for 1y and 10y rate
dev.new()
par(mar=c(2,3,2,.5), mgp=c(2,.6,0), mfcol=c(1,2))
yrange <- data %>% select(-year) %>% range(na.rm=TRUE)
plot(data$year, data$r1, ylim=yrange, col="gray", type="l", lwd=1, xlab="Year", ylab="Percent", xaxs="i", yaxs="i", main="One-year real rate", font.main=1, cex.main=1)
lines(data$year, data$rstar.r1, lwd=3)
lines(data$year, data$rstar.r1.lb, lty=3)
lines(data$year, data$rstar.r1.ub, lty=3)
lines(data$year, data$rstar_mean, lwd=2, lty=5)
legend("topright", c("Real rate", "UC model r*", "Avg. other r*"), col=c("gray", "black", "black"), lwd=c(1,3,2), lty=c(1,1,5))
box()
par(mar=c(2,2.5,2,1))
d2 <- filter(data, !is.na(r10))
plot(d2$year, d2$r10, ylim=yrange, col="gray", type="l", lwd=1, xlab="Year", ylab="", xaxs="i", yaxs="i", main="Ten-year real rate", font.main=1, cex.main=1)
lines(d2$year, d2$rstar.r10, lwd=3)
lines(d2$year, d2$rstar.r10.lb, lty=3)
lines(d2$year, d2$rstar.r10.ub, lty=3)
box()
